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ABSTRACT 



Direct integration leads to the potential gradient of a rectangular plate. Lines with 
equal spacing form a grid for integration over the infinite plane. A source distribution 
is approximated by pattern functions which are centered over each grid point. The 
pattern functions are based upon the sine quotient function. Fourier analysis replaces 
integration in physical space with integration in wave number space. The range of 
integration is limited to lines of finite length in wave number space. The computation 
of potential gradient is provided by subroutines. 
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INTRODUCTION 



The flow around a ship is partly a radial flow from a source distribution and partly 
a circular flow around a vortex distribution. In either case the velocity is the integral 
of the product of an inverse square and the density of sources or vortices. The inverse 
square contributes a spike to the integrands and numerical integration is not feasible. 
A technique for integration through the spike may be found from an analysis of the 
analogous problem in the flow over a plane. 

The potential and the potential gradient of a source or charge which is distributed 
over a plane have applications in hydrodynamics, electromagnetism, and gravitation. 
The field of an irregular three-dimensional source distribution can be constructed 
with source distributions over each of a stack of planes. The potential and the potential 
gradient at a field point near a plane can be computed by two methods. 

In a direct method, the potential at a field point is the integral over the plane of 
the product of inverse distance and the source density, while the potential gradient 
is the integral over the plane of the product of an inverse square of distance and the 
source density. Wherever there is a discontinuity in source density the potential 
gradient is infinite. 

In a Fourier method, the potential is expressed as a solution of Laplace’s equation, 
and the solution is adjusted to fit the boundary conditions. The solution is the sum 
of products of complex exponential functions, and the source distribution is 
approximated as a Fourier series. A Fourier transform replaces integration in physical 
space with integration in wave number space. The time of computation depends upon 
the distance to the field point. There are two methods for the integration in wave 
number space. In the first method, Gauss integration multipliers are applied in a 
roving interval. The time increases with distance. In the second method, integration 
by parts leads to recurrence equations. The time is independent of distance. 



DIRECT INTEGRATION 



Let x,y,z be Cartesian coordinates in a right-handed coordinate system with x and y 
in the plane and with z above the plane. Let cr(x, y) be the source density at a point 
with coordinates x, y. Then the potential is given by the equation 



y, z) = J J 



a(u, v) 

\/(x - u) z + (y - v ) z + 2 2 



dudv 



(1) 



where x,y,z are coordinates of a field point and u, v are coordinates of a source point. 



RECTANGULAR PLATE 



Let a rectangular plate be placed on the plane. Let the length of the plate be 2a, 
and let the breadth of the plate be 2b. Let the origin be at the center of the plate 
with the x-axis in the direction of length and with the y-axis in the direction of 
breadth. 

The potential of a plate with unit source density is given by the equation 



<p{x,y, z) 



n + b-y p+a- 
J — b—y J —a— a 



dudv 

'J u z + v z + z z 



( 2 ) 



Initial integration leads to inverse hyperbolic functions, then final integration is 
completed with integrations by parts. 
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The potential is given by the equation 



, \ . l — i (b-y) , , x . u _, (6 + y) 

<p = (a - x)sinh A 4 (a - x)sinh - 



4 (a 4 xOsinh" 1 
4 (6 - y)sinh -1 
4 (6 + y)sinh _1 
- 2 tan' 



4(a — 


x) 2 + z 2 


(b 


- y) 


V(a 4 


x) 2 + z 2 


(a 


- x) 


\/ (b - 


y) 2 + z 2 


(a 


- x) 



v(6 + y) 2 4 2 2 
(a - x)(6 - y) 



- 2 tan" 



2 \/ (a - x) 2 4(6- y) 2 4 2 2 

t (q + *)(h ~ y) 

2\/( a 4 x) 2 4 (6 - y) 2 4 2 2 



4 (a + x)sinh _1 
4 (6 - y)sinh _1 
4 (6 4 y)sinh _1 
- 2 tan" 1 



4(< 


a - 


x) 2 + z 2 




(6 


+ y) 


4(< 


a 4 


x) 2 + z 2 




(a 


+ x) 


4(6 - 


y) 8 + 2 2 




(a 


+ x) 



\/(6 4 y) 2 4 2 2 

(a - x)(6 4 y) 



2 tan' 



2\/ (a - x) 2 + (6 + y) 2 + 2 2 

( a + x)(6 4 y) 

2 n/(o~T~ c) 2 4(6 4 y) 2 4 2 2 



( 3 ) 



Partial differentiation and cancellation of terms gives the components of the gradient 
of the potential. 

Differentiation with respect to x leads to the equation 



d <P . , _j 

= sinh 

dx 



- sinh" 1 



(t> - y) 



4(a - x) 2 + z 2 

(b ~ y) 



4 sinh 1 
- sinh -1 



( b + y) 



4(a - x) 2 + z 2 

(b + y) 



\/(a + x) 2 + z 2 \/(a + x) 2 + z 2 

differentiation with respect to y leads to the equation 



( 4 ) 



dtp _ 

= sinh 



(a - x) 



dy 



- sinh" 



\/(6 - y) 2 4 2 2 
(a - x) 

V(64 y) 2 4 2 2 



4 sinh' 



- sinh 1 



(a 4 x) 



4(6 - y) 2 + z 2 
(a 4 x) 

\/(6 + y) 2 4 2 2 



(5) 



and differentiation with respect to 2 leads to the equation 



(a-x)(6-y) , A (a — x)(b 4 y) 

— — tan — / • •• — - 4 tan — , - 

& z 2 V( a - x) 2 4 (6-y) 2 4 2 2 2\/(a - x) 2 4 (6 4 y) 2 4 2 2 

(a + x)(6 - y) (a + x)(6 + y) 

+ tan — , ■ ■ ■ + tan — . - - 

zV(a+x) 2 + (6-y) 2 + z 2 zV(a + x) 2 + (6 + y) 2 + z 2 



(6) 



These equations were published in a previous report 1 . They are used in the following 
subroutines. 
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SUBROUTINE RPLTP (AA, AB, AX, AY, AZ, FP) 

********************************************************************************************* 
FORTRAN SUBROUTINE FOR POTENTIAL OF RECTANGULAR PLATE 
********************************************************************************************* 

The half-length a of the plate is given in argument AA, the half-breadth b of the 
plate is given in argument AB, and the Cartesian coordinates x, y, z of the field point 
are given in arguments AX, AY, AZ. The potential of the plate is stored in function FP. 

SUBROUTINE RPLTG (AA, AB, AX, AY, AZ, FX, FY, FZ) 
********************************************************************************************* 
FORTRAN SUBROUTINE FOR POTENTIAL GRADIENT OF RECTANGULAR PLATE 
********************************************************************************************* 

The half-length a of the plate is given in argument AA, the half-breadth 6 of the 
plate is given in argument AB, and the Cartesian coordinates x,y,z of the field point 
are given in arguments AX, AY, AZ. The components of the gradient are stored in functions 
FX, FY, FZ. 



NONUNIFORM PLATE 

A major achievement was the integration when the source density was expressed 
as power polynomials in the surface coordinates. The computation of potential required 
a triplex of recurrence relations. The computation is provided by the subroutine RPLTM. 
The analysis has been presented in a previous report 1 . That RPLTG and RPLTM give the 
same gradients has been verified by computation. 



FOURIER ANALYSIS 



The Fourier transform for a function on a plane is given by the equations 



A(a, (3) = J° j F(x t y) e r ^ ax+ ^v) & y 

F(x, y) = J J A(a, p) e i{ax + 0y) da dp 



( 7 ) 

( 8 ) 



where x, y are Cartesian coordinates in physical space, a,/? are Cartesian coordinates 
in wave number space, F(x t y) is a function in physical space, and A(a, (3) is the Fourier 
amplitude in wave number space. A potential <^>(x, y, z) of a source distribution on the 
plane can be constructed with the equation 



<p(x, y, z) 



JJ 



->P 



A(a, p) e -N“- + /< l»l+««+*v> da d p 



( 9 ) 



which gives a solution of Laplace’s equation wherever z * 0. The vertical component 
of the gradient at z = 0 is given by the equation 



d(p 

dz 



J J Va a + /S a A(a , p) e i(ax+0y) da dp 



( 2 - 0 ) ( 10 ) 



where the sign is 4 - if z > 0 and the sign is - if 2 < 0. In accordance with the Gauss 
theorem the difference in - d<p/dz on opposite sides of the plane is 4rr a(x,y) where 
cr(x,y) is the surface source density at the plane. Thus the Fourier amplitude for an 
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arbitrary distribution of source density is given by the equation 

A(a,p) = f { (j(x,y) e~ i(ax+Pv) dx dy 

2nVcx z + / 3 2 J J 

The range of integration of density is the infinite plane. 
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RECTANGULAR PLATE 



Let a rectangular plate be placed on the plane. Let the length of the plate be 2a, 

and let the breadth of the plate be 26. Let the origin be at the center of the plate 

with the x - axis in the direction of length and with the y - axis in the direction of 

breadth. Let u, v be the coordinates of a source point on the plate and let x,y,z be 

the coordinates of a field point above the plate. 

The potential of a plate with unit source density is given by the equation 



<p(x- y . *) 



■s/JJJ 






2 -/? 2 |*| + i{a(x-u)+/Uy-v)} 



Va 2 + /? 2 

Completion of the integration leads to the equation 



da d(3 du dv 



<p(x. y. z) 



■SJJ 



sin(aa) sin(60) e ~4* z ** z M+H«*+M 



a 



Va 2 + 



da dp 



( 12 ) 



(13) 



What is interesting about this equation is that the distribution of source density over 
a finite plate introduces the factors sin(aa)/a and sin(6/?)//? into the Fourier amplitude. 
The factors concentrate the spectrum of the Fourier integral into bands which are 
centered with respect to the line spectrum of a trigonometric polynomial. 



INTERPOLATION 



The coefficient function for Lagrange interpolation in an infinite set of equally 
spaced data is just the infinite product for the sine quotient function in the equation 



sin u 
u 




* 2 W 



(14) 



It is equal to unity at the origin and is equal to zero wherever else u is a multiple 
of 7T. It is analytic and diminishes with distance from the origin. 

An arbitrary function f(u) is expressed in terms of coefficient functions with centers 
at Am by the equation 



k = -oo 



U - Jc 7T 



(15) 



where f(u k ) is the discrete value of f(u) at the Arth node. 

The product of functions which are centered at kr\ and ran is given by the function 



sin(u - Am) sin(u - ran) 
(u - Am) (u - ran) 



(16) 



Resolution into partial fractions and application of the addition theorem for 
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trigonometric functions leads to the equation 

f +0 ° sin(u-A:7T) sin(u-7n7T) (-1)* +tt T F +0 ° sin 2 (u-A;7T) P + °° sin 2 (u-*m7T) *1 

— 7 — du = — — ; du - du (17) 

J_ TO (u-k 7r) (u-rmr) {k-m)n |_ J _ TO u-kn J _ OQ u-rrnr J 

Inasmuch as the integrands are odd their integrals are zero and the coefficient 
functions are orthogonal with respect to integration. They may be used as the basis 
for approximations by least squares. 

An important property is expressed by the equation 



f 



sin a u 



du = < 



u 



— 7T if a < 0 
0 if a = 0 
l +7T if a > 0 



(18) 



where a is a constant of arbitrary magnitude. Integration by parts leads to the equation 



J + oo O /* + c© . 0 

sin u J sin 2 u 

=— du = du - 7T 

-eo U J-oo U 



(19) 



which gives the normalization factor for least squares approximation. 



SOURCE PATTERN 



A universal source pattern for the infinite plane is given by the equation 



F{x,y) = 




(20) 



The function F(x,y) is unity at the origin and is zero on all equally spaced nodal lines 
which do not pass through the origin. The spacing between nodal lines is 2a in the 
x-direction and is 2b in the y-direction. The function is analytic and diminishes with 
distance from the origin. Application of the Fourier transform and use of the Euler 
theorem introduces the difference of two cosines and introduces the sum of two sines 
into the integrand of the Fourier amplitude. The difference between the cosines is 
zero at the origin and is an even function otherwise. The difference between the cosines 
is eliminated by integration. There remains the sum of sines as expressed by the 
equation 



A(a, /?) 




sin(^-a)x + sin(^+a)x sin(^-/3)y + sin (§+/% 




dx dy 



( 21 ) 



The sines cancel each other during integration wherever a, /S are outside a rectangle 
in wave number space. The interior of the rectangle is located where a,/? satisfy the 
inequalities 



7T 7T 

^ a ^ + — 

2a 2a 



7T 7T 

^ B < + — 

26 H 26 



( 22 ) 



5 



Thus the function F(x,y) is given by the equation 

ab f + 5F r + l; 

F(x, y) = — I I e x(ai+0y) da = ^(x, y) 



and the potential <p(x, y, z) is given by the equation 



tt , n f 

2b r 2a P NT 



, x 2ab f + 2 b r 2 ^ e“ 

<p(x t y, z) = — 

7T J_JL J_JL 

2b 2a 



a 2 +/? 2 |«| + i(ax+/fy) 



Va 2 + /? 2 



da d/? 



(23) 



(24) 



Further integration is possible after Cartesian coordinates a, (3 are replaced by polar 
coordinates /c, 6. The potential is expressed by the equation 



— 



(25) 



Then integration with respect to /c replaces surface integration throughout the interior 
of the rectangle with line integration along the perimeter of the rectangle. 

The three components of the gradient are given by the same formulations except 
for a weighting function w. Thus for the components 



dcp 

dx 

the values of w an 
- i cos 0 

Then each component is the sum 



dcp 

dy 

- i sin 0 
A + B 



dcp 

dz 



+ 1 



(26) 



(2V) 



(28) 



of two integrals where A is the integral along vertical sides of the rectangle and B is 
the integral along horizontal sides of the rectangle. 

For the vertical sides the variables are related in accordance with the equations 



t = tan 6 



(29) 



sin 6 



VTTT* 



cos 0 = 



ViTT* 



a = tc cos 6 = — 

2a 

Let a parameter 6 be defined by the equation 

<5 = ~ j|z|V 1 + f 2 - i(x + yt)\ 

2a 



Then the integral A is given by the equation 



6 f + b 1 - (1 

A = rr - — 

a - 



4 - 6)e 



-6 



w dt 



d6 = + 



dt 



1 4 - t* 



(30) 

(31) 

(32) 

(33) 



for integration along Hoth vertical sides. 
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For the horizontal sides the variables are related in accordance with the equations 

t = cot 9 (34) 



sin 9 = 



1 



VT7t z 



cos 6 = 



t 






(3 = k sin 9 = — 

26 

Let the parameter 6 be defined by the equation 

<5 = ~ {|z|V 1 + < 2 - i{xt + y)t 

cb 

Then the integral B is given by the equation 



B 



a C + a 1 - (1 + 6)e~ 



a f o 

" n b J_b 



w dt 



d8 = - 



dt 



1 + V 



(35) 

(36) 

(37) 

(38) 



for integration along both horizontal sides. 

Wherever |6| S 1 the integrand is replaced by the infinite series in the equation 



-£(* + ■> ( - 4) ‘ 



k = 0 



(Jc + 2)! 



(39) 



Summation of this series is continued until there is no change in the sum. 



GAUSS INTEGRATION 



The integration is completed with the aid of 16-point Gaussian integration 
multipliers. Each interval of integration is subdivided into subintervals and the Gaussian 
integration is applied progressively to successive subintervals. The numbers of 
subintervals for the integration of A and B are the integral parts of the expressions 



2 + 



lyl 

76 




(40) 



Then the Gaussian integration is accurate for any x or y. For large values of x, y the 
time for computation is nearly proportional to the sum of |x| and |y|. 

Formulae for the Taylor series expansion of the source pattern at nodal points in 
a finite grid were used with subroutine RPLTM to obtain gradients at several grid sizes. 
Extrapolation to infinite grid size confirmed the line integration in wave number space. 
The components of the gradient of potential are given by the following subroutine. 

SUBROUTINE RPTNG (AA, AB, AX, AY, AZ, FX, FY, FZ) 

**************************4 1*^********************************^******************************* 

FORTRAN SUBROUTINE FOR POTENTIAL GRADIENT OF RECTANGULAR PATTERN 
**^**^**************************^************************************************************ 



The half-interval a of the pattern is given in argument AA, the half-interval 6 of the 
pattern is given in argument AB, and the Cartesian coordinates x, y, z of the field point 
are given in arguments AX, AY, AZ. The components of the gradient are computed by 
Gauss integration and are stored in the functions FX, FY, FZ. 
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INTEGRATION BY PARTS 



For the integration by parts, the quadratic formula is used to solve the equations 
for 6 to give t as a function of <5. 

For the vertical sides of the rectangle the variable t is given by the equation 



t = 



iyirjT + ix) * \zW + ix) z - (y 2 + z 2 ) 
y z + Z z 



(41) 



Substitution and cancellation lead to the equation 



V i + T z 



\z\(H? + **) ± iy^W + ^) Z ~ ( Z t 

y 2 + z 2 



(42) 



and differentiation leads to the equation 

dt. _ f pyV(^ + **) 2 - (y 2 ± [z[(^ + ix) 

d6 (y z + z z ) [ V(^ + ixf - (y 2 + z 2 ) 



The ± sign in the equations is determined by substitution of the limits of integration 
for t. The sign is + for positive t and is - for negative t . 

Substitution in the equation for A replaces integration with respect to t with 
integration with respect to 6. The natural path of integration in the 2-plane is the 
real axis, whereas the natural path of integration in the (5-plane is an hyperbola with 
vertex where 6 is given by the equation 



(5 = 



7T 

2a 



\\z\ - ix\ 



(44) 



and with axis in the direction of the positive real axis. There is a singularity on each 
side of the imaginary axis. Near the singularity on the positive side (5 is given by the 
equation 



<5 = ~ f Vy 2 + z 2 - ix] + 
2a 



(45) 



Then the limiting values for functions at the singularity are given by the equations 



t = 



iy 






VTT7 2 = 



|Z| 



y + z‘ 



v; 



y + z c 



and by the equation 



V~e- 

dd 






(y z + z z f* 



(46) 



(e - 0) (47) 



The hyperbolic path in the (5-plane is deformed into two straight lines which run from 
the positive singularity to the limits of integration where t = ± a/6. 

For the horizontal sides of the rectangle the variable t is given by the equation 



t = 



ix(^ + iy) ± |z|\/(^ + iy) 2 - (a: 2 + z 2 ) 



x z + z z 



(48) 
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(49) 



Substitution and cancellation lead to the equation 



n/TTT* = 



\z\{^r + iy) ± + iy ) 2 - (z 2 + z 2 ) 



x 2 + z 2 



and differentiation leads to the equation 

dt f r W(^ + iy) 2 -(x 2 +z 2 ) ± I z\(^ + iy) 

d6 (x z + z z ) L + iy ) 2 - (x 2 + Z 2 ) 



(50) 



The ± sign in the equations is determined by substitution of the limits of integration 
for t. The sign is + for positive t and is - for negative t. 

Substitution in the equation for B replaces integration with respect to t with 
integration with respect to 6. The natural path of integration in the f-plane is the 
real axis, whereas the natural path of integration in the 6-plane is an hyperbola with 
vertex where 6 is given by the equation 



6 = 



rr 

26 



IN - iy) 



(51) 



and with axis in the direction of the positive real axis. There is a singularity on each 
side of the imaginary axis. Near the singularity on the positive side 6 is given by the 
equation 






+ 2 2 - iy] + e 



(52) 



Then the limiting values for functions at the singularity are given by the equations 



t = 



xx 



y/: 



Vi + < a = 



X 6 + z c 



y/x* 



X c + 



(53) 



and by the equation 




lz|V! 3 

(x z + Z 2 )f 



(e - 0) (54) 



The hyperbolic path in the 6-plane is deformed into two straight lines which run from 
the positive singularity to the limits of integration where t-± 6/a. 

The weighted derivative 



w 



dt 

dd 



(55) 



is expressed as a power series expansion by 11-point Lagrange interpolation. Let 6 
be given by the equation 

6 = 7} + e (56) 

where 77 is a center of expansion and e is the argument of expansion. Inasmuch as 
the paths of integration are straight lines in the e-plane, all values of e in the range 
of interpolation are given by the equation 

e = ue M (-l^u^+1) (57) 

where e M is the upper limit of interpolation and u is the variable of interpolation. The 
limit e M can be cancelled out of the rational expression for the complex Lagrange 
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coefficient and the derivative can be expanded as a power series in u. In a preliminary 
computation discrete values for u with Chebyshev spacing and discrete coefficients 
were computed for use in a subroutine. Thus the derivatives are expressed as power 
polynomials in the ratio 



u=— (58) 

e M 

while division by e M is deferred until the integration by parts. 

There are three stages of expansion. In the first stage, 77 coincides with the singularity 
and derivatives are expanded as power polynomials in.e 1/2 . The range of interpolation 
is extended into a negative range in order to obtain data for central interpolation. 
The extension in the e 1/2 -plane corresponds to a fold in the e-plane. The interval of 
integration is the upper half of the range of interpolation. In the later stages of 
expansion, 77 is moved further out along the path of integration and the derivatives 
are expanded as power polynomials in e. The interval of integration is the full range 
of interpolation. 

There are tight limitations on the range of interpolation. The radius of expansion 
near the first singularity must be a small fraction of the distance from the first 
singularity to the second singularity. The radius of expansion near the second center 
of expansion must be a small fraction of the distance from the center of expansion 
to either singularity. 

The integration of each series in powers of e 1/2 is completed with a descending 
recurrence. Required for the integration are integrals which are generated by the 
recurrence 



f e^e^de = — e~* + — f 
Jo rn m J 0 



e *e 



1 de 



(59) 



The iteration is started with an initial approximation for 771 = 32. Completion of 
integration is given by integrals which are obtained by the recurrence 




1 - (1 + 6)e~ 6 e m 1 - (1 + 6)e~ s 

7a de = ~ 7 

0 7717} 0 



1 

7717] 




rn - 1 r m 1 - (1 + 6)e 6 _ 

mr) J 0 £ 6 s 



de 



(60) 



where m is alternately half an odd integer and half an even integer. The iteration is 
started with an initial approximation for 771 = 32. Stability of the recurrence is achieved 
by a limitation of |e| to 

The integration of each series in powers of e is cycled in either direction. Required 
for the integration are integrals which are generated by the ascending recurrence 



e m - l e-*de (61) 

The iteration is started with the special integral in the equation 

f e e~ 6 de = e -77 - (1 + e)e~ 6 
Jo 

and the iteration is continued until m = 64. If \e\ ^ 17 then the iteration is cycled in 
descending order with an initial approximation. Completion of integration is given by 



(62) 



J 

Jo 



e m e~ s de = - 



e m e -s 



+ m 



f 

Jo 
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integrals which are generated by the ascending recurrence 



/' 



1 - ( 1 + 6)e~ 



de = 



1 - (1 + 6)e~ 



m - 1 



m 



m : 



e m e 4 de 



_ mr, r m _, 1 ~ (1 

m - 1 Jo 



+ <5)e' 



de 



( 63 ) 



where m is a sequence of integers. The iteration is started with the aid of the special 
integrals in the equations 



r i - (i 

Jo 

C s 1 - (1 + 

Jo u 



+ u)e- u _, , 1-e- 4 " (-6)* +1 

du , 



Jfc = 0 



u)e u 



du = - 1 4- e ~ 6 4 - 74 - log <5 - Ei(-d) = £ 



M) 



*■+ 2 



r *6 <» 

jl - (1 4- u)e~ u )du =- 2 4(5 4 (2 4 6)e' 6 = - £ 
Jo *=0 



*=o (* + 2) 2 fc! 
(-c5) fc+3 



(k 4- 3)(fc + 2)fc! 



(64) 

(65) 

( 66 ) 



and the iteration is continued for 1 1 cycles. If |e| ^ ^ [77 1 the iteration is cycled in 
descending order with an initial approximation. 

In the actual program e is replaced by e/e^, which gives integrals ready for use 
with the Lagrange coefficients, and scales the integrals so that their exponents do not 
exceed the tight limitation on the exponents in IBM computers. The recurrence equations 
can be verified easily by differentiation. 

Although the integration by parts is accurate everywhere that |z| is not small, the 
integration by parts breaks down when y z 4- z 2 is zero along the vertical sides and 
breaks down when x z 4- z 2 is zero along the horizontal sides. Under such circumstances 
the program falls back upon Gauss integration. 

For large values of x,y the time for computation is independent of |x| or \y\. The 
components of the gradient of the potential are given by the following subroutine. 



SUBROUTINE RPTNX (AA, AB, AX, AY, AZ, FX, FY, FZ) 
********************************************************************************************* 

FORTRAN SUBROUTINE FOR POTENTIAL GRADIENT OF RECTANGULAR PATTERN 
********************************************************************************************* 



The half-interval a of the pattern is given in argument AA, the half-interval b of the 
pattern is given in argument AB, and the Cartesian coordinates x, y % z of the field point 
are given in arguments AX, AY, AZ. The components of the gradient are computed with 
integration by parts and are stored in functions FX, FY, FZ. 



GRADIENT OVER GRID 

Exploratory computations have given the gradient of potential on a vertical line 
through each grid point of the source pattern. The gradient has a direction which 
passes through the center and has a magnitude which varies with height. Far from 
the plane the magnitude approaches the full value for the inverse square law, while 
near the plane the magnitude is diminished by shielding by the local structure of the 
source pattern. As distance diminishes toward the center the inverse square of distance 
approaches 00 , while the gradient approaches 2 tt. 
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The potential is given by the equation 



y . z ) 



2a b f + 2b r + z^ e 
7T J _JL J _JL 



+ ^a p'\Z a2+ /* 2 kl+i(ax+/Jy) 



\/a 2 + /3 2 



da d(3 



( 67 ) 



Combinations of the partial derivatives replace the integrand with an exact differential 
which can be integrated at once. Results of integration are given by the equation 



il z 



d£_ x d<p\2ab p" [ _ > /^5|«| +i ( ai+<h/) l + i 
dx dz ) 7T J_JT_ J_J7 L 

2b 2a 



( 68 ) 



The integral is zero if x is a multiple of 2a. Results of integration are given by the 
equation 




The integral is zero if y is a multiple of 2b. Thus the potential gradient has a radial 
direction over grid points where x t y are multiples of 2a, 2b. The integration of the 
inverse square over the surface of the pattern is just the gradient of a point source 
when the field point is over a grid point. 



RECTILINEAR LINE 

When the plane contracts to a line, the potential and the potential gradient become 
symmetric with respect to the line. The position of a field point is located by the 
cylindrical polar coordinates r, 0, z, where r is the distance from the line, 0 is the 
azimuth angle, and z is the zenith distance along the line. The Fourier transform for 
a function on the line is given by the equations 

A(ic) = -^j f(z) e -i ** dz (70) 

J + oo 

Ai^e^dK (71) 

-oo 

where z is the coordinate in physical space and k is the coordinate in wave number 
space. 

A universal source pattern for the infinite line is defined by the equation 



/(z) 




(72) 



where the spacing between nodal points is 2c. Application of the Fourier transform 
and use of the Euler theorem introduces the difference of two cosines and the sum 
of two sines into the integrand of the Fourier amplitude. The difference of cosines is 
an even function of z, and the sum of sines is an odd function of z. Both functions 
are divided by the odd function z. The cosines are cancelled by the integration. There 
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remains the sum of sines as expressed by the equation 



A{ic) 



J_ p sin (^ ~ k)z + sin(§ + k)z ^ 

4 rr J_„ /nz\ 

\ 2c/ 



Then the sines cancel each other whenever /c is outside the range 



( 73 ) 



7T 

2 c 



^ /c ^ + 



2c 



and the amplitude is given by the equation 



A(k) = - 

7T 



(74) 



(75) 



whenever k is inside the range. Thus the function f{z) is given by the equation 



/(z) 



7 T 




2 c 




cos(/cz) dtc 



( 76 ) 



where the range of integration in wave number space is finite. 

The Laplacian in cylindrical polar coordinates is given by the equation 



V-Vcp 



d z w 1 dcp 1 d z <p d z cp 

_| — -f — -| — 

dr 2 r dr r z d<f) Z dz z 



(77) 



The solutions of Laplace’s equation are given in many texts. In this case the appropriate 
solution contains the modified Bessel function of the second kind, because this function 
diminishes to zero with increasing argument 2,3 
The potential is given by the equation 



2c r* 

c p{r t z) = — e XKZ K 0 {fcr) d/c 

7T J_jl 

2c 



_ 4C 

” Jo 



cos (/cz) K 0 (/cr) die 



(78) 



where K 0 {/cr) is the Bessel function of order zero. The Bessel functions satisfy the 
equation 



d 

— K 0 (xr) = - k K x (fer) (79) 

dr 

where K^/cr) is the Bessel function of order one. The Bessel functions satisfy the limits 

K 0 (fcr) -► - log(/cr) K t {ter) -♦ — (r ^ 0) (80) 

ter 

at small radius. Thus the potential satisfies the limit equation 

dw 

- 27t r — =47 t f(z) (r 0) (81) 

dr 

where f(z) is the linear density a(z) in accordance with the Gauss theorem. 
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components 


of the 


gradient of potential 


are 


given by the equations 






2c 


r* . , x 


4c 


IT 

Czc 




v 'f' 


4 - 


ic e tKZ K 1 (/cr) die = + 




ic cos(icz)K i (ter) die 


(82) 


dr 


7T j 


Ljl 

2c 




'o 


d(fi 


2 ic 


r** , . 


4c 


7 r 

Czc 








ic e Kn(icr) d/c - + 




ie sin (iez)K 0 (ier) die 


(83) 


dz 


7T , 


2c 


7T , 


Jo 



Inasmuch as these integrals are limited to a finite range, they are evaluated with 
Gauss integration multipliers. 

The potential gradient of a rectilinear source pattern is given by the following 
subroutines. 

SUBROUTINE BSSLK (MO, AZ, IN, FK) 

****+***+**********************************+************************************************* 
FORTRAN SUBROUTINE FOR MODIFIED BESSEL FUNCTION OF INTEGRAL ORDER 
*************++*************+***********************+^**+********^*****************++******** 

The mode of operation is given in MO. The real and imaginary parts of the argument z 
are given in array AZ, and the integer order n is given in IN. The modified Bessel 
function of the second kind is computed by series expansion, rational approximation, 
and recurrence relations. If MO = 0, the real and imaginary parts of the function K n (z) 
are stored in array FK. If MO = 1, the real and imaginary parts of the function e z K Tl {z) 
are stored in array FK. 

SUBROUTINE LPTNG (AC, AR, AZ, FR, FZ) 

FORTRAN SUBROUTINE FOR POTENTIAL GRADIENT OF LINEAR PATTERN 

^^^*^^^H:^^^*^t******************************************************************************* 

The half-interval c of the pattern is given in argument AC, and the cylindrical 
coordinates r, z of the field point are given in arguments AR, AZ. The components of 
the gradient are computed by Gauss integration and are stored in the functions FR, FZ. 

DISCUSSION 



The ubiquitous function 



sin x 
x 



(84) 



is sometimes called the diffraction function because it appears in the diffraction of 
waves at an aperture, but it appears in other applications as well. It should have a 
name of its own. It is proposed that it be called the sine quotient function. 

The source density is a scalar while the vortex density is a vector. The sine quotient 
function can be used for the interpolation of the source density or either of the 
components of vortex density. Then the potential gradient of the pattern function is 
used directly with source density in the computation of radial flow, or in a vector 
product with either component of vortex density in the computation of circular flow. 

A technique for the computation of ideal flow over arbitrary bodies has been invented 
by Hess and Smith 4,5 . The same technique has been adopted for computation of flow 
over ships by Dawson and Dean 6,7 . In this technique the surface of the body is subdivided 
by a grid, then each panel of the grid is projected onto a plane which spans the four 
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corners of the panel. It is assumed that the quadrilateral projection has a uniform 
source density. Each plane quadrilateral is assembled from plane triangles, for which 
the potential gradient is given by analytical expressions. 

In the case of an infinite plane, uniform quadrilaterals would be reduced to 
rectangular plates. The potential gradient of a square plate is compared with the 
potential gradient of a source pattern in the following table. 



Potential Gradients 

a - b = 0.5 x = y = 0 



z 


Inverse 

Square 


Plate 


Pattern 


0 


oo 


2tt 


2n 


0.0625 


256 


5.58063622 


5.41516931 


0.125 


64 


4.90438059 


4.68192879 


0.25 


16 


3.70918087 


3.53400628 


0.5 


4 


2.09439510 


2.09491828 


1.0 


1 


0.80543168 


0.86196571 


2.0 


0.25 


0.23543002 


0.24793473 


3.0 


0.11111111 


0.10812127 


0.11106356 


4.0 


0.0625 


0.06154089 


0.06249869 


5.0 


0.04 


0.03960461 


0.03999996 


6.0 


0.02777777 


0.02758643 


0.02777777 


7.0 


0.02040816 


0.02030466 


0.02040816 


8.0 


0.015625 


0.01556424 


0 01562500 



The entries in the table were computed with the aid of RPLTG, RPTNG, and RPTNX. The 
table shows that the inverse square law gives accurate gradients for z > 4. 

For rectangular plates the field point must be kept away from the edges, where 
the gradient is infinite, whereas for pattern functions there is no restriction on the 
location of the field point. 



CONCLUSION 

Pattern functions which are based upon the sine quotient function are useful for 
the computation of the potential and the potential gradient of source distributions 
in the plane. Directly over a grid point the potential gradient has a direction which 
passes through the center of the pattern, but has a magnitude which is shielded from 
the full inverse square law. 
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